!!             mod_2dflu.f90
!! This is the module containing the defination and commons for the 
!! serial code. As all the serial ffts requal arrays of same type 
!! we expect this should not be changed as we go from one machine  
!! to another.
module mod_2dflu
  use omp_lib
  use technical
  implicit none
  save
  public
  !! ----velocites  -------------
  double precision,target,allocatable,dimension(:,:) :: psi,omega, & 
       jac_old,ukx,uky 
  double precision,pointer :: ux(:,:),uy(:,:)
  double precision,allocatable,dimension(:,:) :: omega_tm
  double precision, allocatable, dimension(:) :: xx,yy
  double precision :: xx0=0.,yy0=0.,xmax=2*pi,ymax=2*pi
  !! ---parameters of simulations -----------------
  double precision,allocatable,dimension(:,:) :: E_Omega
  integer :: pmax
!! ------logical variables -----------
  logical :: lfluid=.true.,lparticle=.false.
  logical :: limmersed_boundary=.false.
  logical :: lforcing=.true.,lfourier_forcing=.true.,lrspace_forcing=.true.
  logical :: loutputpy=.true.
  logical :: lwrite_velocity_output=.true.,lascii_output=.true.,lwrite_vel_rspace=.true.
  logical :: lwrite_streamfunction=.false.
  character (len=labellen) :: datadir='data/'
!! -----common parameters ----------------------------
  integer :: thousand,nn,maxiter,nrun,n1,n2,n1h,n1hf, &
       nshell,navg,ksqr_max,mm,iinner,nloop,nalias
  double precision :: tmax
  integer :: kini,kdrag,kasqr,no_of_fmode,t_n,ip
  character(len=labellen) :: hydro_ini='zero'
  double precision :: vis,vis2,delta,length,factor,mu,mu2,finp,oamp=0.01
  double precision :: fixed_energy1,fixed_energy2,fixed_energy3
  logical :: lcalc_Gradu=.false.
  real*8::dx,dy
!! -----global arrays --------------
  double precision,allocatable,dimension(:) :: time_increment,alpha,diss_rate,tot_energy,&
    energy_time,enstrphy_dissrate,energy_dissrate
  integer,allocatable,dimension(:) :: den_state
  integer,allocatable,dimension(:,:) :: fmode 
!! ---- For fourier transform -----------------
  integer,dimension(2) :: dim
  double precision :: scale
!! ---energy injection --------------------------
  real*8::fw,fu
  double precision,allocatable,dimension(:,:) :: loc1
!! ---parameters of simulations -----------------
  double precision :: edr,tlr_micro_scale,dissipation_scale,vkrms,Rlambda, &
    tau_eddy,integral_scale,Rbox,Rintscale
!! -----------For FFTW --------------------------------------
!! ---------------------------------------------------------
  real*8,allocatable,dimension(:,:)::dx_ux,dx_uy,dy_ux,lambda
!! ---------------------------------------------------------
contains
!! *********************************************************** !
subroutine intialize_grid()
  integer :: ix,iy
  write(*,*) 'Initializing grid'
  n1 = nn !the whole boxsize
  n2 = nn 
  n1h = n1/2
  n1hf = n1/2 + 1
  ksqr_max = 2*(nn/2)*(nn/2) 
  nshell = int(1.414*nn/2.0d0) + 1 	
  nalias = int(nn/3) 
  kasqr = 2*nalias*nalias
!	write(*,*) kasqr,nn
!! nshell = \sqrt{ 2 nn}, the diagonal of the box in fourier space. 
!	write(*,*) nshell
!! The length of the box is fixed to be twice \pi.  
  length = 2.0d0*pi
  factor = 2.0d0*pi/length  
  dx = length/dble(n1);
  dy = dx
  allocate(xx(n1))
  allocate(yy(n2))
  do ix=0,n1-1
    xx(ix+1)=xx0+dx*ix
  enddo
  do iy=0,n2-1
    yy(iy+1)=yy0+dy*iy
  enddo
  open(unit=3,file='grid.temp',status='unknown')
  do ix=1,n1
    write(3,*) ix,xx(ix)
  enddo
  do iy=1,n2
    write(3,*) iy,yy(iy)
  enddo
  call write_grid()
!
endsubroutine intialize_grid 
!! ---------------------------------------------------------
subroutine write_grid()
  if (loutputpy) then
    open(unit=11,file='grid_data.py',status='unknown')
    write(11,9) 'def rgrid_param():'
    if (lfluid) then 
      write(11,12) 'lfluid=1'
    else
      write(11,12) 'lfluid=0'
    endif
    if (lforcing) then 
      write(11,12) 'lforcing=1'
    else
      write(11,12) 'lforcing=0'
    endif
    if (lparticle) then 
      write(11,12) 'lparticle=1'
    else
      write(11,12) 'lparticle=0'
    endif
    if (limmersed_boundary) then 
      write(11,12) 'limmersed_boundary=1'
    else
      write(11,12) 'limmersed_boundary=0'
    endif
    write(11,12) '# ------------------------------- #'
  else
    open(unit=11,file='grid_data.ascii',status='unknown')
  endif
  write(11,10) 'n1=',n1
  write(11,10) 'n2=',n2
  write(11,11) 'length=',length
  write(11,11) 'xx0=',xx0
  write(11,11) 'yy0=',yy0
  write(11,11) 'xmax=',xmax
  write(11,11) 'ymax=',ymax
  write(11,11) 'dx=',dx
  write(11,11) 'dy=',dy
! 
! This last line must remain consistent with the rgrid in the file
! grid.py
!
  if (loutputpy) & 
    write(11,12) 'return lfluid,lforcing,lparticle,limmersed_boundary,n1,n2,length,xx0,yy0'
  close(11)
9 FORMAT(A)
10 FORMAT(T4AI8)
11 FORMAT(T4Af18.10)
12 FORMAT(T4A)
endsubroutine write_grid
!! ---------------------------------------------------------
subroutine init_uu(omega_amp,kinput)
  integer, intent(in) :: kinput
  double precision, intent(in) :: omega_amp
  integer :: i1,i2
  select case(hydro_ini)
    case('none','zero')
      omega=0.
      psi=0.
    case('kolmogorov-x')
      call set_kolmogorov_x(omega_amp,kinput)
    case('random-omega')
      call set_random_omega(omega_amp,kinput)
    case default 
      call fatal_error('init_uu','velocity initial condition not chosen')
    endselect
endsubroutine init_uu
!! ---------------------------------------------------------
subroutine set_kolmogorov_x(omega_amp,kinput)
  integer, intent(in) :: kinput
  double precision, intent(in) :: omega_amp
  integer :: ix,iy
  double precision :: kf,x
  write(*,*) ' Kolmogorov flow as initial condition, x dependent'
  write(*,*) 'omega = oamp*kin*sin(kin*x)'
  write(*,*) 'with kin=',kinput,'  oamp=',omega_amp
  kf = dfloat(kinput)
  do iy=1,n2
      omega(1:n1,iy) = omega_amp*kf*sin( kf*xx(:) )
  enddo
endsubroutine set_kolmogorov_x
!! ---------------------------------------------------------
subroutine set_random_omega(omega_amp,kinput)
  double precision,dimension(3) :: ran
  double precision :: rk2,rk,ek1,vk,p11,p12,p31,p22,p33,p23,rk2inv
  double precision :: iniamp,omegak,n_sqr,rkini,mabs,prim,pert,x,y
  integer :: i1,i2,i3,k1,k2,k3,ksqr,mshl,ik,iseed,ireal,iimag,m1,m2
  integer,parameter :: nsize=1
  integer ::ns
  integer,dimension(nsize) :: seed
  integer,dimension(3) :: tarray
  integer,dimension(2) :: tofday
  integer :: ierr
  double precision :: rphi1,rphi2, komega
  double precision, intent(in) :: omega_amp
  integer, intent(in) :: kinput
!! -----------------------------------------------------------------
!! vorticity has a flat spectra with random phase 
  write(*,*) ' Rnadom omega as initial condition: '
  write(*,*) 'with kin=',kinput,'  oamp=',omega_amp
  komega=dfloat(kinput)
  n_sqr = dfloat(n1*n2)
  call itime(tarray)
  seed(1) = tarray(3)*9765 + tarray(2)*1734 + tarray(3)
!  write(*,*) 'iseed-',iseed
  call init_random_seed()
  call random_number(rphi1); call random_number(rphi2)
  rphi1 = rphi1*2*pi - pi
  rphi2 = rphi2*2*pi - pi
!! --------- 
  do i2 = 1,n2
    y = yy(i2)
    do i1 = 1,n1
      x = xx(i1)
     !
      prim = dcos(komega*y+rphi1) + dcos(komega*x+rphi2)
      pert = 0.0d0
      do m2 = 0,2!Nx/2
        do m1 = 0,2!Ny/2
          mabs = dsqrt(dble(m1*m1 + m2*m2))
          !
          if(m1**2+m2**2.ne.0)then
            pert = pert + &
              1.0d-4*(dsin(m1*x+m2*y) + dcos(m1*x+m2*y))*dble(m2*m2)/mabs
          endif
          !
        enddo
      enddo
!        omega(i1,i2) = (prim + pert)*vis
      omega(i1,i2) = (prim + pert)*omega_amp*komega
    enddo
  enddo
!--------------------------------------------
endsubroutine set_random_omega
!! ---------------------------------------------------------
subroutine init_random_seed 
  INTEGER :: i, n, clock
  INTEGER, DIMENSION(:), ALLOCATABLE :: seed
  n=1
!  
  CALL RANDOM_SEED(size = n)
  ALLOCATE(seed(n))
!  
  CALL SYSTEM_CLOCK(COUNT=clock)
 ! 
  seed = clock + 37 * (/ (i - 1, i = 1, n) /)
  CALL RANDOM_SEED(PUT = seed)
!  
  DEALLOCATE(seed)
END subroutine init_random_seed
!! *********************************************************
end module mod_2dflu
